Solubilities and Transport Properties of CO2, Oxalic Acid, and Formic Acid in Mixed Solvents Composed of Deep Eutectic Solvents, Methanol, and Propylene Carbonate

Recently, deep eutectic solvents (DES) have been considered as possible electrolytes for the electrochemical reduction of CO2 to value-added products such as formic and oxalic acids. The applicability of pure DES as electrolytes is hindered by high viscosities. Mixtures of DES with organic solvents can be a promising way of designing superior electrolytes by exploiting the advantages of each solvent type. In this study, densities, viscosities, diffusivities, and ionic conductivities of mixed solvents comprising DES (i.e., reline and ethaline), methanol, and propylene carbonate were computed using molecular simulations. To provide a quantitative assessment of the affinity and mass transport of CO2 and oxalic and formic acids in the mixed solvents, the solubilities and self-diffusivities of these solutes were also computed. Our results show that the addition of DES to the organic solvents enhances the solubilities of oxalic and formic acids, while the solubility of CO2 in the ethaline-containing mixtures are in the same order of magnitude with the respective pure organic components. A monotonic increase in the densities and viscosities of the mixed solvents is observed as the mole fraction of DES in the mixture increases, with the exception of the density of ethaline-propylene carbonate which shows the opposite behavior due to the high viscosity of the pure organic component. The self-diffusivities of all species in the mixtures significantly decrease as the mole fraction of DES approaches unity. Similarly, the self-diffusivities of the dissolved CO2 and the oxalic and formic acids also decrease by at least 1 order of magnitude as the composition of the mixture shifts from the pure organic component to pure DES. The computed ionic conductivities of all mixed solvents show a maximum value for mole fractions of DES in the range from 0.2 to 0.6 and decrease as more DES is added to the mixtures. Since for most mixtures studied here no prior experimental measurements exist, our findings can serve as a first data set based on which further investigation of DES-containing electrolyte solutions can be performed for the electrochemical reduction of CO2 to useful chemicals.


INTRODUCTION
During the past few decades, carbon capture, utilization, and storage (CCUS) technologies have been in the spotlight of the academic and industrial research as a means for reducing the concentration of CO 2 in the atmosphere. 1 One promising CCUS route is the utilization (e.g., reduction) of CO 2 as a feedstock for the production of value-added products. 2,3 Several technologies are available for the reduction of CO 2 , e.g., photocatalytic, thermal, and electrochemical. Electrochemical processes have distinct advantages such as the lack of complex reaction pathways, cost-efficiency, and relatively high reduction efficiencies. 4−6 CO 2 can be electrochemically converted to a number of valuable materials and fuels, spanning polymers, acids, alcohols, and hydrocarbons. 3,7 Valuable CO 2 electroreduction products include formic and oxalic acid, which are the simplest forms of monocarboxylic and dicarboxylic acids, respectively. 8,9 The CO 2 electoreduction to these acids require only two moles of electrons per mole of product and have a high market price. 7,10 Formic acid can be produced with high Faraday efficiencies (>95%) and current densities (150 mA cm −2 ) using gas diffusion electrodes. 10 In 2018, formic acid was reported to have a total market value of $756.5 MM, with a market price of approximately $400/tonne. Formic acid is mostly used in agriculture, the production of leather and textiles, and in the pharmaceutical industry. 11 Oxalic acid is mainly used in the pharmaceutical and textile industry. 12,13 Oxalic acid has a global market value of $715 MM and a market price of ca. $500/tonne. 14 Despite the tremendous progress that has been made in the field of electrochemical processes during the past few decades, significant challenges and limitations still remain. 6 The main challenges are the high overpotential requirements and the low selectivity toward the desired products. To overcome these limitations, many factors have to be considered when designing and optimizing an electrochemical conversion process, e.g., the electrochemical cell configuration, catalyst, and type of electrolyte. 7,15 The role of the electrolyte is of particular importance since it constitutes the medium for the conversion reactions and controls the transport of the different chemical species to the catalysts. 15 Consequently, selecting the optimum electrolyte/solvent for a conversion process can enhance the performance of electrochemical conversion processes. 16 To this purpose, many electrolytes have been tested through the years, e.g., aqueous and organic solvents and ionic liquids (ILs). 5 ILs have been considered for these processes due to high thermal stability, ionic conductivity, and absorption of CO 2 . The use of ILs has also been shown to reduce the required overpotential and undesirable side reactions in electrochemical conversions, while the ILs themselves can act as a co-catalyst. 17,18 Deep Eutectic Solvents (DES) are an emerging class of solvents sharing similar properties and advantages with ILs. 19−28 Many DES, e.g., choline-based, can be easily prepared from mixing naturally occurring substances and, thus, are cheaper to produce than most ILs. 29,30 Compared to ILs, the use of DES in electrochemical processes is not as widely investigated. High viscosities can be a limiting factor toward application of pure DES as electrolytes for the electrochemical reduction of CO 2 . 31 To exploit the benefits of DES in such processes while overcoming the drawbacks, mixing DES with other solvents has been considered. Vasilyev et al. 31 showed that the CO 2 reduction reaction takes place in the presence of various choline-based DES, such as reline and ethaline (which are formed by mixing choline chloride with urea and ethylene glycol, respectively, in the ratio 1:2). Vasilyev et al. 31 also observed that the efficiency of CO 2 reduction increased upon the addition of DES in the originally used electrolyte, i.e., ethylene glycol.
A first approach for examining the feasibility of solvents containing DES in electrochemical applications is to investigate the thermo-physical properties of these solvents and of the respective mixtures with the reactants and products. For example, the solubility and diffusivity of solutes (e.g., CO 2 , products) in electrolytes are very important properties since they often are limiting factors in electrochemical conversions. Excess properties and solubilities of the solutes in the solvents are equally important for, e.g., the design of downstream separation processes following the conversion of CO 2 to the value-added products. While experiments are traditionally used to measure properties of fluid mixtures, molecular simulations are less costly and, therefore, can assist in the initial screening of a large number of solvents for electrochemical processes.
Molecular simulation also provides the necessary fundamental understanding of the physical/chemical mechanisms at the atomistic scale. For these reasons, molecular simulations have been widely used to compute various properties relevant to electrochemical applications. 32−40 In this work, the solubilities and self-diffusivities of CO 2 , oxalic acid, and formic acid in mixtures of DES with organic solvents are computed by means of Monte Carlo (MC) and Molecular Dynamics (MD) simulations. Self-diffusivities, densities, viscosities, and ionic conductivities of the solvent mixtures are also computed as a function of the composition of the mixtures. Two DES, i.e., reline and ethaline, are considered here. The organic solvents considered are methanol and proplyene carbonate. These solvents have been used as electrolytes for the conversion of CO 2 to formic acid and oxalic acid, respectively. [41][42][43][44] Our study shows that the reline− methanol mixtures have slightly lower affinity toward CO 2 and that the addition of DES to the organic solvents increase the solubilities of oxalic and formic acids. The densities and viscosities increase with the mole fraction of DES, except for the density of ethylene-propylene carbonate (due to the higher density of the pure organic component compared to the DES). In contrast, the self-diffusivities of all molecular species vastly decrease due to the increasing viscosity. For all mixed solvents, the ionic conductivities show a nonmonotonic behavior with the DES content. Initially, the ionic conductivity increases until a maximum value, and then a sharp decrease is observed as more DES is added. This behavior is in line with prior studies on aqueous DES solutions, and reline−ethaline mixtures. 26,45,46 Overall, comparisons of our simulation data with the limited available experimental measurements are in reasonable agreement. This paper is organized as follows. In Section 2, the computational details regarding the MC and MD simulations and the force fields used are provided. The results of the thermodynamic and transport properties are presented in Section 3. In the same section, an analysis of the hydrogen bonding behavior of the system is performed. The conclusions of this study are discussed in Section 4.

METHODS
Molecular simulations are performed for the following solvents: methanol, propylene carbonate, reline, ethaline, and mixtures of ethaline-propylene carbonate, ethaline-methanol, and reline-methanol. The mole fraction of DES in the different mixtures is defined as follows: HBD HBA HBD HBA organic (1) where N HBD , N HBA , and N organic is the number of hydrogen bond donors, acceptors, and organic molecules. For example, in the case of ethaline-methanol mixtures, N HBD , N HBA , and N organic correspond to the total number of ethylene glycol, choline chloride, and methanol molecules, respectively. 2.1. Force Fields. Nonpolarizable force fields consisting of bonded (i.e., bond streching, angle bending, and torsions) and nonbonded (i.e., Lennard-Jones and Coulombic) terms were used to simulate all species in this work. The TraPPE force field was used to model CO 2 47 and methanol. 48 For oxalic acid, the modified OPLS force field proposed by Doherty and coworkers 49,50 was used. Formic acid was modeled using the modified OPLS force field parametrized by Salas et al. 49,51 which yields improved predictions for the dielectric constant. Propylene carbonate parameters were taken from the work of Silva and Freitas, 52 who adopted GAFF and refitted the charges. The DES were modeled using the GAFF 53 force field consistently with our previous studies. 21,24,26,27,54 For choline, urea, and ethylene glycol, 1−4 interactions were scaled by a factor of 0.5 for both Lennard-Jones and Coulombic interactions. The charges of choline chloride were scaled by a factor of 0.8 and 0.9 in reline and ethaline, respectively. 55,56 This implementation yields accurate predictions for various thermophysical properties of DES as shown by Perkins et al., 55,56 Salehi et al., 24 and Celebi et al. 26,27,57 The Lennard-Jones interaction parameters between unlike species were computed using the Lorentz−Berthelot combining rules. 58 All force field parameters and the functional forms of the bonded and nonbonded terms used in this study are available in the Supporting Information.
2.2. Monte Carlo Simulations. In this work, MC simulations were performed to compute the excess chemical potentials (μ ex ) and Henry coefficients (H), which are used to quantify the solubilities of solutes (i.e., CO 2 , oxalic acid, and formic acid) in different mixed solvents. For a component i, the excess chemical potential μ i ex follows from 59 where μ i and μ i IG are the chemical potentials of the component and the ideal gas at the same conditions, respectively. For a specific solute−solvent combination, μ i ex indicates the affinity of the solute toward the solvent as it is related to the activity coefficient γ i of component i: 60,61 γ ρ where ⟨ρ 0i ⟩ is the ensemble average number density of pure component i, ⟨ρ i ⟩ is the ensemble average number density of i, x i is the mole fraction of i, μ 0i ex is the excess chemical potential of pure i with respect to the ideal gas, k B is the Boltzmann constant, and T is the temperature in units of K. The Henry coefficient of the solute, H i is defined as 59 where P i and f i are the partial pressure and fugacity of the solute, respectively. H i is directly related to μ i ex as follows: 62 where ρ is the number density of the solvent and T is the temperature in units of K. All MC simulations were carried out with the open-source software package Brick-CFCMC, 63,64 which utilizes the Continuous Fractional Component (CFC) method 65,66 (i.e., gradual insertion/deletion of fractional molecules during the simulations). The degree of interaction between a fractional with the surrounding molecules is varied using a scaling parameter λ (0≤ λ ≤ 1), which is a degree of freedom in an expanded ensemble formulation. 67 For more details on the CFCMC method, the reader is referred elsewhere. 65−69 Recently, a thermodynamic integration feature has been developed in Brick-CFCMC for computing μ ex based on 64 where U is the energy of the system, and the brackets ⟨···⟩ denote an ensemble average. During CFCMC simulations, separate scaling parameters are used for intermolecular Lennard-Jones and electrostatic interactions. The scaling parameters are continuous functions of λ and are implemented such that electrostatic interactions are not switched on before fully scaling down the Lennard-Jones interactions. For more details, including the scaling functions, the reader is referred to the work of Polat et al. 64 To compute μ ex of CO 2 , oxalic acid, and formic acid in different solvents using eq 6, the λ space was discretized into 50 bins. Separate simulations in the NPT ensemble were performed for each solute with a fixed value of λ to compute λ ∂ ∂ U . Subsequently, numerical integration of eq 6 was performed. More details on the thermodynamic integration feature of Brick-CFCMC can be found in the recent work of Polat et al. 64 μ ex and H were computed for mixtures with 0 ≤ x DES ≤ 0.4 at 298.15 K and 1 atm and for pure reline and ethaline at 350.15 K and 1 atm. A cutoff radius of 12 Å was used for both the Lennard-Jones and the Coulombic potential in all MC simulations except for the ones of pure DES in which a cutoff radius of 10 Å was used. Electrostatic interactions were handled with the Ewald summation method with a relative precision of 10 -6 . During the MC simulations, trial moves were selected with the following probabilities: 35% translations, 35% rotations, 29% changes in the internal configuration of molecules (i.e., angles and dihedrals), and 1% volume changes. A minimum of 8 × 10 5 cycles were carried out for equilibration and 8 × 10 5 cycles for production. At each MC cycle, the number of the trial moves performed equals the number of molecules of the system.

Molecular Dynamics
Simulations. MD simulations were performed for the computation of the densities, number of hydrogen bonds (HBs), shear viscosities, and self-diffusion coefficients. All MD simulations were carried out using the large-scale atomic/molecular massively parallel simulator (LAMMPS). 70 The initial configurations were generated with the PACKMOL package. 71 Long-range electrostatic interactions were handled using the particle−particle particle-mesh (PPPM) method with a relative error of 10 −6 . The cutoff radius was set to 12 Å for both Lennard-Jones and the shortrange part of the Coulombic interactions. Periodic boundary conditions were imposed in all directions. The Verlet algorithm with a time step of 1 fs was used to integrate Newton's equations of motion. Temperature and pressure were maintained constant using the Nose−Hoover thermostat and barostat with coupling constants of 100 and 1000 fs, respectively.
Transport properties were computed with the OCTP (onthe-fly computation of transport properties) plugin in LAMMPS 72 which yields the mean-squared displacements (MSDs) of dynamical properties as a function of time. The transport coefficients can be then obtained by linear regression to the long-time MSDs at time-scales where the slopes as a function of time are equal to 1 in a log−log plot. Diffusion coefficients are computed from 58,72 The Journal of Physical Chemistry B pubs.acs.org/JPCB Article where D i MD is the self-diffusivity of species i, r j,i (t) is the position vector of the j th molecule of species i at time t, and N i is the number of molecules of species i in the system. The shear viscosity η follows from 58,72 (9) where V is the volume of the system, P αβ ′ are the components of the traceless pressure tensor, P αβ are the off-diagonal components of the pressure tensor, and δ αβ is the Kronecker delta. All self-diffusion coefficients were corrected for finite-size effects using the Yeh-Hummer (YH) equation: 74−76 where D i is the corrected self-diffusion coefficient corresponding to the thermodynamic limit, η is computed from MD simulations and does not depend on the system size, 77 where e is the elementary charge and q i is the charge of the molecules of species i. Eq 11 has been shown to be a relatively good approximation for obtaining the ionic conductivities of ionic species in a computationally efficient way. 26,54,80,81 For all mixtures considered here, only the charges and the selfdiffusivities of choline and chloride were used in the NE equation since the rest of the species are charge-neutral (i.e., HBDs and the organic solvents). The ionic conductivity can also be computed using the appropriate Green−Kubo and Einstein relations (i.e., cross correlation of charge fluxes/ displacements). 79 The MD simulations of the solvents with x DES ranging from 0 to 1 were performed at 298.15 K and 1 atm. A list of the solvents studied here and the number of molecules used for each species is shown in Table 1. For the computation of the self-diffusivities of CO 2 , oxalic acid, and formic acid in the different solvents, five solute molecules were used. This helps to drastically improve the sampling of MSDs, while it practically corresponds to infinite dilution. The MD simulation scheme was as follows. Initially, an energy minimization using the conjugate-gradient method with a tolerance of 10 −4 was performed. Then, equilibration runs in the NPT ensemble were carried out for 10−20 ns, depending on x DES . Finally, production runs of 10−100 ns were carried out in the NVT ensemble from which all properties were computed. For each system, averages and standard deviations were computed over 5 independent MD simulations, each one starting from a different initial configuration. Visual molecular dynamics (VMD) 82 was used for the HB analysis. The criterion for the formation of a HB was a cutoff distance of 3.5 Å between the donor and acceptor atoms and a cutoff angle of 30°between the donor-hydrogen-acceptor atoms. 83    The Journal of Physical Chemistry B pubs.acs.org/JPCB Article all systems with maximum absolute deviations of 1.2%, 1.1%, and 1.0% for reline-methanol, ethaline-methanol, and ethalinepropylene carbonate mixtures, respectively. These low deviations serve as validation of the accuracy of the selected force fields.
As expected, the densities of the methanol-containing solvents increase considerably with the addition of DES, due to the large difference between the densities of the pure components (i.e., the densities of methanol, ethaline, and reline are 778, 1120, and 1200 kg/m 3 , respectively). Relinemethanol mixtures are denser than ethaline-methanol mixtures for any x DES . This is also expected since the density of pure reline is higher than that of pure ethaline. Ethaline-propylene carbonate mixtures have higher densities compared to the methanol-containing ones for x DES ≤ 0.8. In these systems, the density decreases with the addition of DES (opposite behavior from the methanol mixtures); however, this decrease is not large. The density of ethaline−propylene carbonate mixtures decrease by 5% as x DES increases from 0 to 0.8. This is mainly due to the similar densities of pure ethaline and pure propylene carbonate. As shown in Figure 1, no experimental data are available for the ethaline-propylene carbonate mixtures for x DES > 0.2. On the basis of the excellent agreement between the MD and experiments for x DES < 0.2 and for the rest of the ethalinecontaining solvents, our new predictions can be considered trustworthy.

Excess Chemical Potentials and Henry Coefficients.
In this section, we present the computed excess chemical potentials and Henry coefficients of CO 2 , oxalic acid, and formic acid in different solvents consisting of a DES (i.e., reline or ethaline) and an organic cosolvent (methanol or propylene carbonate). Our approach was verified by comparing the solubility computed from MC simulations with experimental measurements for the case of CO 2 in pure methanol at T = 313.15 K and 2 atm. Figure 2 shows the values of the average derivative of the energy with respect to the λ parameter in the CFCMC simulations as a function of λ. Using thermodynamic integration (eq 6), we obtain μ CO 2 ex = −3.27 kJ/mol, and from this we compute H CO 2 = 0.58 MPa. This value deviates by around 4% from the respective experimental Henry coefficient reported by Xia et al. 85 This small deviation indicates that the chosen force fields and the method are reliable.
The computed values for μ ex and Henry coefficient of CO 2 , oxalic acid, and formic acid in the different solvents are listed in Table 2. As can be seen, the solubility of CO 2 in pure methanol and pure propylene carbonate is almost equal (absolute deviation of ca. 5%). Clearly, the addition of DES in these organic solvents reduces the CO 2 solubilities. For x DES = 0.4, the solubilities of CO 2 are reduced by ca. 50% and 30% for reline-methanol and ethaline-methanol mixed solvents, respectively. For the same x DES in ethaline-propylene carbonate mixture, the solubility of CO 2 is reduced by ca. 20%. The Henry coefficients listed in Table 2 indicate that solvents containing ethaline are slightly better adsorbents of CO 2 than reline-containing mixtures.
Interestingly, the Henry coefficient of oxalic acid in DESmethanol mixtures is much lower compared to the one in the pure solvents (i.e., x DES = 0 and 1). We speculate that this could be due to an interplay between hydrogen bonding interactions and a commensurate fit of the oxalic acid molecule in the liquid structure. Overall, the computed Henry coefficients show that adding a choline-based DES to the organic solvent increases the solubilities of oxalic acid and formic acid. While the solubility of CO 2 is reduced as a result of adding a DES to methanol or propylene carbonate, it is important to note that the reduction is not very large and the Henry coefficients are still relatively high. The mixed solvents investigated here have higher CO 2 Henry coefficients compared to aqueous solutions at the same conditions, which are typically used in electrochemical processes. 5 This is an important finding since the design of an electrolyte with high CO 2 solubility could potentially improve conversion rates by increasing the concentration of CO 2 at the surface of the electrode. 15 3.2. Transport Properties. 3.2.1. Viscosities. The computed viscosities of the ethaline-propylene carbonate, ethaline-methanol, and reline-methanol mixtures are shown in Figure 3 as a function of the DES content. Available experimental data 86−88 are also shown in this figure along with the available experiments. Clearly, the viscosity increases with the DES content. Interestingly, for the ethaline-propylene carbonate mixture, this is the opposite behavior compared to the densities discussed earlier. The reason is that the density of pure propylene carbonate is slightly higher compared to pure ethaline, while the viscosity of propylene carbonate is much smaller than that of ethaline. This can be mainly attributed to the fact that, unlike pure propylene carbonate, pure ethaline has a strong hydrogen bonding network (this is discussed in detail in the following section). The viscosities of the pure organic solvents are predicted with deviations from experiments of ca. 7% and 20% for methanol and propylene carbonate, respectively. For all mixed solvents, the deviation between the computed and the experimentally measured viscosities increases with the addition of DES in the mixtures (absolute standard deviations range from ca. 7% to 44% and ca. 20% to 50% for ethaline-methanol and ethaline-propylene carbonate, respectively). Although these deviations seem rather high, the available experimental data are limited (e.g., no experimental data exist for the reline-methanol mixture) and the uncertainties in the computed values are quite large, ranging from 7 to 25% (due to the difficulty in sampling the slow dynamics caused by the relatively low temperature). It is also important to note that large deviations are reported  25 In absolute values, the predicted viscosities from MD simulations are satisfactory, while the qualitative behavior of the systems is captured accurately. Given the scarcity of experimentally measured viscosities for most of the mixtures considered here, our MD data can serve as a first set of predictions to aid the design of industrial processes and further motivate experimental efforts. To improve the accuracy of the computations, further modifications to the force fields, combining rules, and/or charge scaling should be considered. Such an investigation is beyond the scope of the present study. As can be seen in Figure 3, reline-methanol mixtures are significantly more viscous than ethaline-methanol for the whole DES composition range. For x DES = 0.8, the viscosity of reline-methanol is higher than the viscosity of ethalinemethanol by almost a factor of 3. This is not surprising since pure reline is significantly more viscous than ethaline (i.e., η reline = 455 mPa s and η ethaline = 62 mPa s at 298 K). For x DES < 0.2, the viscosities of both mixtures are within the same order of magnitude. In the range (x DES ≤ 0.6), ethaline-propylene carbonate mixtures exhibit the largest viscosity. For x DES > 0.8, reline-methanol viscosities are the highest, and ethalinepropylene carbonate viscosities become comparable to the viscosities of ethaline-methanol.
Overall, our results reveal that the addition of DES to the organic solvents have a very strong effect on the viscosities. From a practical point of view for electrochemical processes, this finding dictates the careful selection of the composition of the mixtures since large viscosities can limit mass transport and, thus, reduce the current density of the electrolytes. 31 To this end, DES with relatively low viscosities such as reline or ethaline (or other) can be promising.
3.2.2. Self-Diffusivities. As slow diffusion rates can be a limiting factor in electrochemical processes, it is essential to integrate an electrolyte that yields sufficient mass transfer of the reactants and products to and from the catalyst. 4,15,89 Since no experimental diffusivity data are available for the mixtures studied here, our results are the first step toward the screening of solvents for an optimum electrolyte containing DES, methanol, and propylene carbonate. In this section, we present the computed self-diffusivities of all the species in the DES− organic solvent mixtures, and the self-diffusivities of infinitely diluted solutes (CO 2 , oxalic acid, and formic acid) in these mixtures.
The computed self-diffusion coefficients of the different molecular species in the reline-methanol, ethaline-methanol, and ethaline-propylene carbonate mixed solvents are shown in Figure 4. All reported diffusivities were corrected for finite-size effects using eq 10. As can be clearly seen, the self-diffusivities of all components monotonically decrease as the DES composition increases. This is mainly due to the increasing viscosities of the mixtures upon the addition of DES (see Figure 3), resulting in reduced mobilities of the different species. Due to the very high viscosity of propylene carbonate, the ethaline-propylene carbonate mixtures are the most viscous for x DES < 0.8. This is clearly reflected to the self-diffusivities of all species in this mixture for the same range of DES compositions, which have lower values compared to the  The molecular weight (MW) and the hydrodynamic radius are known factors to greatly affect the diffusivity of a molecule in a solvent. 90 In DES (and DES-containing mixtures), the presence of an extended network of HBs is another crucial factor affecting mass transport. 26,56,57 Choline, which is the heaviest species (MW ≈ 104.2 g/mol) among all HBDs and HBAs, has the lowest diffusion coefficient in all mixtures and DES compositions. Interestingly, the diffusivity of the much lighter chloride (MW ≈ 35.5 g/mol) is comparable to that of choline and lower than the diffusivities of both the HBD (i.e., urea and ethylene glycol with MW of 60.06 and 62.07 g/mol, respectively). In ethaline-containing mixtures, the diffusivity of ethylene glycol is higher than that of chloride by ca. 28%. In reline-methanol, the diffusivity of urea is slightly higher than that of chloride. Similar trends for the diffusivities of the HBD and HBA species were observed in the study by Celebi et al. 26 for aqueous DES mixtures. This behavior can be explained by the HB network within the DES. As suggested by Perkins et al., 56 the fact that urea diffuses faster than most of the components in reline (despite having almost twice the MW of chloride) can be attributed to the formation of many HBs with other urea molecules and the anions. This can be clearly seen in Figure 5a, in which the computed HBs between the components of the DES are shown. Due to the varying number of molecules used in the MD simulations of different solvents (see Table 1), the number of HBs were normalized to represent a system containing 100 DES molecules. The number of the organic molecules follows from x DES . As shown in Figure 5, in the methanol-containing solvents all HB combinations monotonically increase as more DES is added to the mixture. The number of HBs formed between the various species increases ca. 2 to 6 times in the range of x DES = 0.1−1.
In the reline-methanol mixture, the rise in the number of urea−urea HBs is impressive, going from 25 to 101 (per 100 reline molecules). In the ethaline-methanol mixture, the anion-HBD HBs are also significantly increased, ranging from 14 to 86 (per 100 ethaline molecules) in the range of x DES = 0.1−1.
In the same mixture, the HBs between the HBD molcules are more than quadrupled (10 to 42/100 DES). The gradual development of this strong HB network is the main reason for the increasing viscosities and decreasing diffusivities of the different species in the methanol-containing solvents discussed earlier. In ethaline-propylene carbonate, the numbers of HBs formed between the various species do not significantly vary The computed number of HBs formed between the organic solvents and the DES species are shown in Figure 6. Again, the number of HBs is normalized to represent a system containing 100 methanol or propylene carbonate molecules. The number of HBD and HBA follows from x DES . In contrast to propylene carbonate, methanol can form HBs with all the DES components (and with other methanol molecules). Thus, as x DES increases, the methanol-methanol HBs are being depleted, and methanol forms HBs with the HBDs and HBAs. As can be clearly seen in Figure 6a,b, methanol primarily forms HBs with the HBD (urea or ethylene glycol) and secondarily with the anions. This HB behavior, combined with the relatively low MW of methanol (≈32 g/mol), are the main reasons for the fast self-diffusivities shown in Figure 4d. The lack of HBs between propylene carbonate and most of the DES components can be seen directly in Figure 6c and indirectly in Figure 5c. In the latter, the absence of competition between the organic component and the DES species to form HBs is the main reason for the almost constant HBs numbers between the HBA and HBD of the ethaline, with the only exception being the increasing HBD−HBD HBs. The self-diffusion coefficients of infinitely diluted CO 2 , oxalic acid, and formic acid in the different solvents are shown in Figure 7 as a function of x DES . Consistently with our findings for the solvents, the diffusivities of all solutes decrease as the DES mole fraction increases. In all mixtures, CO 2 has the highest self-diffusivity followed by formic acid and oxalic acid. This order is in line with the molecular weights of these solutes. The highest diffusivities of all solutes are observed in ethaline-methanol. For x DES < 0.4, all solutes diffuse faster in the methanol-based solvents than in the ethaline-propylene carbonate mixture. As discussed earlier, this can be mainly attributed to the high viscosity of the ethaline-propylene carbonate mixture for this composition range. At x DES = 0.6, the lines representing the self-diffusivities of CO 2 (Figure 7a   The Journal of Physical Chemistry B pubs.acs.org/JPCB Article slowest due to the fact that this mixture is the most viscous one in this concentration range as shown in Figure 3. Because of the very small number of solutes used in the MD simulations (corresponding to infinite dilution), a solute-DES or soluteorganic solvent HB analysis is not a very accurate descriptor for explaining the diffusivity behavior of the solutes, thus these HBs are not reported here.

Ionic Conductivities.
Another important property to optimize when designing electrolytes for electrochemical applications is ionic conductivity since electrolytes ensuring fast electron transfer are essential for high-performance electrochemical conversions. Recently, ionic liquid-based electrolyte solutions have been studied for the electroreduction of CO 2 to valued-added products. 5 To the best of our knowledge, no experimental data are available for the ionic conductivities of the DES-organic solvent mixtures considered here. The ionic conductivity of pure reline has been measured experimentally by various groups to be in the range of 0.024− 0.764 S m −1 for T = 293−353 K, respectively. 45,46,91,92 It is important to note that the actual values reported in literature significantly vary depending on the experimental technique used and the purity of the DES. For example, Agieienko and Buchner 45 reported an electric conductivity of 0.024 S m −1 for pure reline at 298 K, while at the same conditions, Mjalli and Ahmed 46 report a value of 0.18 S m −1 , which is an order of magnitude higher. Celebi et al. 26 reported a value of 0.09 S m −1 computed in MD simulations at 303 K. Here, a value of 0.11 S m −1 has been computed for T = 298 K. The measured ionic conductivity of ethaline ranges from ca. 0.62 to 2.08 S m −1 in the temperature range of 293−353 K. 46,92 At room temperature it is equal to ca. 0.70 S m −1 (the exact value depends on the experimental study). Here, we computed a value of 0.63 S m −1 , which is in reasonable agreement with the experiments. Since a thorough validation of the computed conductivities for the mixtures of DES with methanol and propylene carbonate is not possible due to the absence of experimental measurements and due to the fact that the NE equation has been shown to slightly overpredict conductivities, 54,80,93 our results should be interpreted mostly qualitatively.
The computed ionic conductivities of all mixed solvents studied in this work are shown in Figure 8 as a function of x DES . For all solvents, the ionic conductivities exhibit a nonmonotonic behavior. As x DES increases, the ionic conductivities initially increase until a maximum value after which a sharp decline is observed. This can be explained by the fact that as ionic content (i.e., DES) is added to the mixture, the ionic conductivity initially increases up to the maximum value. However, the sharp increase of the viscosity due to the formation of the strong HB network within the DES (see Figures 5 and 6) restricts the mobility of the ions, causing the decline of κ after a certain x DES . This nonmonotonic behavior is fully consistent with the MD data by Celebi et al. 26 and the experiments by Agieienko et al. 45 for aqueous reline and ethaline solutions. Mjalli and Ahmed 46 also observed this behavior for reline-ethaline mixtures.
Methanol-containing solvents have higher ionic conductivities compared to ethaline-propylene carbonate. For x DES ≤ 0.6, this difference is significant, i.e., a factor of 2 to 6. The only exception is for x DES = 0.8, for which the ethaline-propylene carbonate solvent exhibits slightly higher ionic conductivity than the reline-methanol one. This is in-line with the viscosity of these mixtures, which follows the exact same trend. The    Figures 4 and 7), the peak of ionic conductivity for ethalinepropylene carbonate is shifted toward higher x DES . As x DES approaches 1, the hydrogen-bonding network in the DES becomes extensive, causing the viscosity to significantly increase and, thus, the ionic conductivities of all solvents to reach their minimum. The only exception is ethaline-propylene carbonate, due to the very high viscosity of the pure organic component.

CONCLUSIONS
The electrochemical reduction of CO 2 to value-added products, such as formic and oxalic acid, is considered to be a promising carbon utilization route for partially mitigating the greenhouse effect. Recently, DES have been considered as possible electrolytes for the reduction reactions of CO 2 as a nontoxic and cost-efficient alternative to ionic liquids. Despite the distinct advantages of these solvents, the applicability of pure DES as electrolytes is hindered by high viscosities.
Mixtures of DES with organic solvents can be a promising way of designing superior electrolytes by exploiting the advantages of each solvent type. In this study, the Henry coefficients and self-diffusivities of CO 2 , oxalic acid, and formic acid in relinemethanol, ethaline-methanol, and ethaline-propylene carbonate mixed solvents were computed using MC and MD simulations. The densities, viscosities, self-diffusivities, and ionic conductivities of the mixed solvents were also computed. The simulations were performed at T = 298 K, P = 1 atm, and mixture compositions x DES = [0,1]. Our simulations showed that the Henry coefficients of CO 2 in the ethaline-methanol and ethaline-propylene carbonate mixtures are in the same order of magnitude as the pure organic components. The reline-methanol mixtures were found to have slightly lower affinity toward CO 2 . Overall, the addition of DES to the organic solvents was found to increase the solubilities of oxalic and formic acids. The densities and viscosities of the mixed solvents monotonically increase with the mole fraction of DES. The only exception was observed for the density of ethalinepropylene carbonate which shows the opposite behavior due to the fact that the pure organic component is much denser than the pure DES. The self-diffusivities of all components (i.e., HBDs, HBAs, methanol, and propylene carbonate) in the mixtures significantly decrease as the mole fraction of DES increases. Interestingly, the self-diffusivities of the infinitely diluted CO 2 and oxalic and formic acids decrease by 1 to 2 orders of magnitude as the composition of the mixture shifts from the pure organic component to pure DES. Our HB analysis revealed that the number of HBs between the DES species is vastly affected by the presence of methanol. As the mole fraction of DES increases, the HBs formed between methanol molecules are being depleted and methanol starts forming new HBs with the HBAs and HBDs of reline. At the same time, a sharp increase in the HBD-HBD and HBD-anion is observed. In sharp contrast, the presence of propylene carbonate has a smaller effect on the HB network of the DES, since it cannot form HBs with most of the DES species. A nonmonotonic behavior was observed for the computed ionic conductivities as a function of composition, which initially increased with the mole fraction of DES, showed a peak at a specific mole fraction for each mixture, and then decreased as more DES was added to the mixture. This finding is in-line with prior literature studies on aqueous DES solutions and mixtures of reline with ethaline. From an application point of view, the thermophysical data produced in this study suggests that the mixtures with low DES content could be the most practical in electrochemical processes since these mixtures exhibit lower viscosities compared to pure DES, higher ionic conductivities than the pure organic solvents, and good absorption capabilities. For most of the mixtures studied here, no prior experimental measurements exist, thus our findings can be considered a first approach based on which further experimental and theoretical studies of DES containing electrolyte solutions can be performed.
Force field parameters and tabulated raw data for the computed thermo-physical properties (PDF)